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SUMMARY 

Minimum altitude-loss flight of a sailplane through a given vertical wind 
distribution is discussed. The problem is posed as an optimal control problem, 
and several numerical solutions are obtained for a sinusoidal wind distribution. 


INTRODUCTION 

The problem of determining the optimal sailplane trajectory through a pre- 
scribed vertical wind distribution for minimum altitude loss is formulated and 
solved as an optimal control problem. The flight is assumed to take place in 
a vertical plane over a fixed range, and the rotational or pitch dynamics of 
the sailplane are neglected. Sailplane lift coefficient serves as the control 
function in the nonlinear point-mass equations of motion. 

For oscillatory vertical wind distributions, this problem belongs to the 
class of "optimal dolphin soaring" problems. In qualitative terms, these 
problems exhibit solutions for which the sailplane speed is decreased in 
upcurrents to prolong the altitude gain and increased in downcurrents to lessen 
the altitude loss (ref. 1). Earlier solutions to these problems have assumed 
either piecewise-static flight (equilibrium glide through segments of constant 
vertical wind — see, for example, reference 2) or quasi-static flight (kinematic 
equations of motion only — see, for example, references 3 through 6). Thus, 
the primary distinguishing feature of this paper is the use of the full nonlinear 
translational equations of motion and the corresponding use of a modern optimal 
control algorithm for numerical solutions. Additional research on the appli- 
cation of optimal control theory to dynamic sailplane performance problems may 
be found in references 7 through 9. 


PROBLEM FORMULATION 

A brief derivation of the equations of motion used here is provided in 
the Appendix. The basic assumptions are: flight in a vertical plane, uniform 

gravity acceleration g and atmospheric density p, a point-mass sailplane of 
constant mass m, and vertical wind of magnitude W. If the vertical wind dis- 
tribution is further assumed to be independent of altitude (W = W(x)), then 
the right-hand sides of the equations of motion do not depend on altitude Y. 
The altitude equation (lOA) can therefore be Incorporated into the performance 
index (altitude loss) 
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J = Y(0) - Y(tj) 


i 


(- Y)dt 


- (w 


+ V sin y)dt 


and need not be regarded as a differential constraint. 


( 1 ) 


Furthermore, it will be convenient to regard the range X as the independent 
variable rather than the time t. Since the final range, X(tj) = X^, is to be 

specified, this change of variables will result in a fixed "end-time" optimal 
control problem which is inherently easier to solve than a variable "end-time" 
problem. The range equation (9A) can also be omitted from consideration as a 
differential constraint. It must be tacitly assumed, however, that the optimal 
trajectory will not Include any kind of looping maneuver which would result in 
zero values for V cos y. Using (9A) then, the performance index (1) becomes 


J = 



W + V sin Y 
V cos y 


dX 


and the remaining equations of motion, (llA) and (12A), become 


( 2 ) 


dX 


- I pV CpS/(2m) + [(V cos y) (dW/dX) + glsin y | /(V cos y) 
^ “ I pVCj^S/(2m) - [cos y(dW/dX) + g/V]cos y | /(V cos y) 


respectively. 


(3a) 

(3b) 


Finally, the nondimensional quantities 


X - X/Xj, V = V(gXj)"^, w = W(gXj)“^ (4) 

are introduced. The resulting optimal control problem may be stated as follows. 
Find that control function u(x) , 0 £ x £ 1, which minimizes the augnented 
performance index 




V sin y + w(x) 
V cos Y 

1 


dx + 


K-l/ [V/ 


''stall ■ 


'‘2' L ' 


1 - v/v ] ^dx 
max 


subject to the second-order dynamic system 


^ - [nCp(u)v^ + (1 + w)sin y]/(v cos y) , v(0) - v 


(5) 
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(6a) 



(6b) 


^ = [nC^(u)v^ - (1 + w)cos Y]/(v^ cos y). Y(0) = Y^ 

and subject to the terminal state constraints 


- v(l)/v -1 = 0 
■L 0 

(7a) 

1(^2 = Y(1) - Yq = 0 

(7b) 

2 

where C^(u) = + a^ 

(8) 

C^(u) = (2 sin^ u - 1) 

(9) 

max 


n = jP (S/mg)gX^ 

(10) 


and where w(x) is the prescribed wind distribution, w = (dw/dx) (dx/dt) 

= (dw/dx) V cos Y, and is the fixed range. Minimum altitude-loss equilibrium 

glide (still air) values are adopted for the fixed and equal initial and terminal 
state values, and Several additional explanatory comments are required. 

First, note that minimum (stall) and maximum (flutter) state Inequality 
constraints on the airspeed are enforced using integral interior penalty 
functions (ref. 10) shown in terms two and three, respectively, of equation 
(5). Thus, a sequence of optimal control problems (3) - (10) must be solved 
for specified positive penalty constants and K 2 . The penalty constants are 
then increased between subproblems. The solution obtained from each subproblem 
is used as starting data for the subsequent subproblem. The sequence of sub- 
problems is terminated when each penalty function value is sufficiently small. 

Secondly, it may be observed that the lift coefficient is bounded via the 
transformation (9). That is, for any value of the control funcuion u(x) , the 
control inequality constraints 

- < C^Cu) < ( 11 ) 

max max 

are satisfied where is a specified constant. Also, note that a quadratic 

drag polar (8) is used?®*More accurate drag polars may be used instead provided 
that an analytical relation is available between C and Cj^. Finally, it should 
be observed that for a fixed wing loading, the nondimensional aerodynamic 
parameter n in equation (10) is proportional to the specified terminal range X^. 

This is the basic timal control problem considered here. A variation of 
this problem will also b presented later. 
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NUMERICAL RESULTS 


All computations have been performed on coupled IBM 360/65 and Itel AS/5 
computers using a FORTRAN IV compiler and double precision arithmetic. The 
numerical integration of the required differential equations has been performed 
using a standard fourth-order Runge-Kutta method with 100 fixed uniform 
integration steps. 

The numerical results have been obtained for the case of a sinusoidal wind 
distribution 


w(x) = w^ sin(2iTx), 0 _< x ^ 1 (12) 

and the Nimbus II open-class sailplane using the gradient projection algorithm 
presented in reference 11. The sinusoidal wind distribution (12) is simply 
an idealized model of an oscillatory vertical wind which satisfies a "continuity" 
condition: the integral of w(x) over the fixed range is zero. The values 

for the coefficients 


a^ = 0.009278, 


32 = - 0.009652, 


a^ = 0.022288 


(13) 


of the quadratic drag polar (8) for the Nimbus II have been obtained from a 

least-squares fit of data taken from the manufacturer's velocity polar. For 

standard sea level conditions and a wing loading, mg/S, of (32) (9. 81) N/m^, 

the aerodynamic parameter in equation (10) is given by n = 0.01916 X^. 

Additional constant data chosen include: C, = 1.4, (gX^)^ v ^ ,, =18 m/s 

L t stall 

and (gX,)^ v =70 m/s. Finally, in qualiPt^ive terms, the gradient projec- 
tion method Ts a direct method in the sense that the control function u(x) 
is changed during each iteration so as to produce both a decrease in the per- 
formance index value (eq. (5)) and full satisfaction of the terminal state 
constraints (eqs.(7)). 


Specified Initial State 

In this case, the initial and final state are to be held fixed and equal. 

In particular, the values 

(gX.)^ V = 28.1676 m/s and y = “ 0.019106 rad (14) 

to o 

are to be used in equations (6) and (7) and correspond to the minimum altitude- 
loss equilibrium glide conditions for the Nimbus II in still air with a drag 
polar given by equations (8) and (13) . The vertical wind amplitude is chosen 
as (gX )^ w. = 2 m/s, and a fixed range of 1000 meters is used. 

L Pi 

The resulting optimal trajectory and the corresponding optimal lift 
coefficient distribution are presented in figures 1 and 2, respectively. The 
optimal flight can be divided into three successive segments: an initial climb, 

a maximum C arc, and a dive followed by a short pull-up. The initial climb 

L* 

l9 intuitively reasonable since the sailplane must gain as much altitude as 
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possible while In the Initial upcurrcnt. The maximum arc Is a continuation 
of the first phase and lasts as long as the wind Is strong enough to sustain 
It. The following dive Is made to pass through the downcurrent as quickly as 
possible. The final pull-up is necessary to meet the terminal state constraints 
(7) . The stall speed inequality constraint was active for this solution, but 
the maximum speed constraint was not. 

The minimum altitude loss for this optimal trajectory is only 12.19 m. 

By comparison, the minimum altitude loss during an equilibrium glide in still 
air over the same 1000 m range is 19.11 m. This represents a 36% altitude - 
loss reduction. 


Free but Equal Initial and Final States 

Here, the initial and final speed and local flight path angle values are 
no longer specified, but the respective initial and final values are still 
required to be equal. The gradient projection algorithm, as described in 
reference 11, can accomodate the addition of the two control parameters v(0) 
and y(0) representing variable initial states. However, the presence of these 
same two control parameters in the terminal state constraints necessitates a 
further modification to the projection operator equations. 

Since the optimal trajectories are now being selected from a larger class, 
additional performance gains are expected. However, for a final range of 1000 
m, the minimum altitude loss improves only 1.4%: from 12.19 m to 12.01 m. 

For comparison purposes, the optimal trajectory is also shown in figure 1. 

In this case, the optimal trajectory exhibits a higher altitude gain during 
the climb phase, a longer maximum arc and a lesser altitude loss in the 
downcurrent when compared with the previous solution. The initial (and final) 
airspeed has increased approximately 1.2 m/s to 29.384 m/s. 


Effects of Wind Amplitude 

The initial and figal states are again free but equal. For an increased 
wind amplitude of (gX^)'- ^ m/s, a substantial improvement is obtained as 

may be noted from the optimal trajectory shown in figure 3. A net altitude 
gain of 5.158 m is now available over the 1000 m course. Clearly, for higher 
amplitudes of an oscillating vertical wind, more energy can be extracted from 
the wind to sustain cross-country flight. 


Effects of Varying the Fi.xed Range 

The wind amplitude is held fixed a- 5 m/s, and f ce but equal initial and 
final states are again considered. Changes in tlie final range X^ affect only 
the constant aerodynamic parameter n and the characteristic length and time 
used in the nondimensionallzation. Varying X^ is equivalent to varying the 
frequency of the sinusoidal wind distribution for sustained flights. 
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Upon reducing from 1000 m to 500 m, a radically different optimal tra- 
jectory was obtained. The corresponding optimal trajectory and optimal 
distribution are shown in figures A and 5, respectively, and will be referred 
to as a Type II solution. In this case, a dive-first climb-later flight pattern 
is observed. The average speed is much higher, and the net altitude gain Is 
considerably higher than for the earlier Type 1 trajectory. The net altitude 
gain of 23.10 m exceeds that achieved on the previous Type I solution for double 
the range. For the Type II trajectories, the maximum speed Inequality constraint 
is active rather than the stall speed constraint. 

By employing the results of previous Type I solutions as starting data, 
it was possible to also obtain a Type I solution for = 500 m. However, this 
Type I relative minimum solution yields a net altitude loss of 4.45 m. Thus, 
at least for relatively small X and large W , the Type II solution is decidedly 
superior to Type I solutions. As a matter or conjecture, the Type II solution 
may cease to exist for sufficiently large final ranges and/or for sufficiently 
small wind amplitudes. Two additional solutions were also obtained: a Type 

I solution for = 750 m and a Type II solution for = 625 m. 


Effect of Wing Loading 

If the nominal wing loading, mg/S, is increased by 15%, the aerodynamic 
parameter 0 becomes 0.01666 X^ (see equation (10)). Nothing else is changed. 

In this case, a Type I solution was obtained for X^ = 1000 ra, = 5 m/s and free 
but equal initial and final states. The resulting^optimal trajectory provides 
a net altitude gain of 1.14 m which is nearly 4 m less than the comparable 
5.16 ra obtained earlier for the nominal wing loading. Both the optimal 
trajectories and the optimal lift coefficient histories for the two solutions 
are very similar. 

The key results for the eight optimal solutions presented here are 
summarized it table I. 


CONCLUSIONS AND DISCUSSION 

The optimal control problem treated here is of at least moderate difficulty 
in view of the state variable inequality constraints present. Relatively few 
numerical solutions are currently available because of the considerable compu- 
tational effort involved. However, several tentative conclusions emerge from 
the computational results obtained thus far. 

1) For a sinusoidal vertical wind distribution, which serves as a simple 
model of a zero range-averaged oscillatory wind, substantial altitude 
savings are available when compared with optimal equilibrium glides 
in still air. The relative advantage increases for higher wind 
amplitudes . 
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2) Equal Initial and final state vector elements can be treated as addi- 
tional control parameters in the optimal control algorithm and there- 
fore varied as part of the optimization process. The additional 
altitude gains obtained in this case are, however, rather small. 

3) For relatively short ranges and high wind amplitude, it is possible 

to obtain optimal trajectories which exhibit an unexpected "dive first, 
climb later" maneuver sequence. Optimal trajectories of this second 
type involve higher speeds and better final altitude gains than the 
usual "dolphin" style optimal trajectory. 

4) As expected, an increase in sailplane wing loading increases the min- 
imum altitude loss when other conditions are held fixed. 

The most surprising finding of this study is the apparent existence of 
two distinct types of extremal solutions at least for restricted ranges of the 
parameters involved. Clearly, the Type II trajectory deserves further research 
effort . 

There is also a need to obtain results for other wind distributions and 
to make definitive comparisons with previous optimization studies which do not 
incorporate the full translational equations of motion. The related problem 
of minimum-time flight through a given vertical wind distribution for a specified 
altitude loss is of perhaps even greater interest. Research on this latter 
problem is currently underway. 
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APPENDIX: DERIVATION OF THE EQUATIONS OF MOTION 

The sailplane velocity vector, relative to the surrounding air, is 
j given by (see figure 6) 

V - R - Tf (lA) 

vrtiere R is the Inertial velocity vector for the sailplane, and - [0, W(X,Y)] 
is the deterministic vertical wind distribution. The angle y shown in figure 
6 is the usual flight path angle for the case of no wind. In the inertial 
(X,Y)-coordlnate frame, the translational equations of motion for unpowered 
flight in a uniform gravity field may be written as 

m^ = T+D+G = m(V + W ) (2A) 

where L and D are the usual aerodynamic lift and drag forces, respectively, G 
is the weight force, and m is the constant sailplane mass. Equation (2A) can 
be rearranged to yield 

• « 

V = (L + ^ + G) /m - W 

= L/m + D/m + [0,-g]^ [0, -X(3W/3X) - YOW/3Y)]^ (3A) 

where g is the constant gravity acceleration. 


Primarily because of the definitions of lift and drag, it is desirable 
to rewrite equation (3A) with respect to a rotating (5,n)-coordinate frame 
defined by the unit vectors e., directed toward V, and e , directed normal 
to V along the lift vector. TOe inertial time derivativS of V, referred to 
this rotating frame, is then given by 


• 

V = (dV/dt) , + w X V 

(C.n) 

= V e^ + (y e^) ^ (V e^) 

= V e^ + V y e^ (4A) 


where e = e^ 
in the? rotating 


e . Note that any vector in the plane, say A, can be expressed 
(^»R)-frame using the following rotation matrix. 


A 




A = I A^ j = I cosy siny 

-siny cosy 

Using equation '5A), one obtains 

G/m - vT = - [XOW/3X) + Y(3W/3Y)+ g][siny e^. + cosy e^] 

The term, g + W » g + X(3W/3X) + Y(3W/3Y), Is often referred to as 
gravity acceleration. 


(5A) 

(6A) 

the apparent 
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Finally, using equations (4A) and (6A) and the usual expressions for 
lift and drag, equation (3A) can be written in scalar form as follows. 

V - - pV^CjjS/(2m) - [XOW/3X) + Y(3W/3Y) + g]sln y (7A) 

Vy - pV“Cj^S/(2m) - [X(3W/3X) + Y(3W/3Y) + g]cos Y (8A) 

From equation (lA) and the inverse of equation (5A), Che kinematic equations 
may be written as 

X =■ V cos Y 

Y = W + V sin Y UOA) 

Therefore, equations (7A) and (8A) beco e 

Y = - pV^CpS/(2m) - [(V cos y)OW/3X) 

+ (W + V sin y)( 3W/3Y) + g]sln y (UA) 

Y =• pVCj^S/(2m) - [cos yOW/ 3X) + (W/V + sin y)OW/3Y) 

+ g/V]cos Y (12A) 
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fixed boundary conditions 
wing loading Increased 13% 


ALTITUDE 

Y(0) - V(l) CHANGE 

rad. m 


-0.0191® 

-12.187 

0.0511 

-12.012 

-0.0528 

5.158 

-0.6286 

23.098 

-0.5285 

11.283 

-0.0089 

-4.454 

0.0039 

-4.452 

0.0903 

1.14o' 
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Figure 1.- Optimal trajectories of Type I: * 2 m/s, 

Xf » 1000 m, Vstall ■ 18 m/s, CL ma-, » 1.4. 



RANGE, X, m 

Figure 2.- Optimal lift coefficient distribution for fixed and equal boundary 
conditions: • 2 m/s, Xf ■ 1000 m, Vstall " 1® m/®t ^I^nax " 
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Figure 3.- Optimal Type I trajectory for high wind amplitude: 
Xf « 1000 m, " 18 
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I Figure 4.- Optimal Type II trajectory: - 5 m/a, - 500 a 
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■ 5 m/s. 


V--* • 70 m/s. 
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